
% PSO算法 
function [x_best, y_best, z_best] = pso(fitness_func, N, iter_max)
c1 = 2.0; % 学习因子常数
c2 = 2.0;
w_max = 0.9; % 不变性权重上限
w_min = 0.4; % 不变性权重下限

% 定义搜索范围
x_min = -4;
x_max = 4;
y_min = -4;
y_max = 4;

% 初始化粒子位置、速度、适应度及最佳位置
x = x_min + rand(1, N)*(x_max-x_min);
y = y_min + rand(1, N)*(y_max-y_min);
vx = randn(1, N);
vy = randn(1, N);
z = zeros(1, N);
for i = 1:N
    z(i) = fitness_func(x(i), y(i));
end
z_best = z;
x_best = x;
y_best = y;

% PSO算法迭代
for iter = 1:iter_max
    w = w_max - (w_max-w_min)*iter/iter_max; % 不变性权重
    for i = 1:N
        % 更新速度
        vx(i) = w*vx(i) + c1*rand()*(x_best(i)-x(i)) + c2*rand()*(x_best(randi(N))-x(i));
        vy(i) = w*vy(i) + c1*rand()*(y_best(i)-y(i)) + c2*rand()*(y_best(randi(N))-y(i));
        
        % 更新位置
        x(i) = x(i) + vx(i);
        y(i) = y(i) + vy(i);
        x(i) = min(max(x(i),x_min),x_max);
        y(i) = min(max(y(i),y_min),y_max);
        
        % 计算适应度
        z_new = fitness_func(x(i), y(i));
        
        % 更新最佳位置
        if z_new < z(i)
            z(i) = z_new;
            x_best(i) = x(i);
            y_best(i) = y(i);
            if z_new < z_best
                z_best = z_new;
            end
        end
    end
end

% 输出最佳位置及最小值
fprintf('The optimal solution is (%.4f, %.4f), and the minimum value is %.4f.\n', x_best(find(z==min(z))), y_best(find(z==min(z))), min(z));

% 绘制函数三维图
[X,Y] = meshgrid(x_min:0.1:x_max,y_min:0.1:y_max);
Z = fitness_func(X,Y);
surf(X,Y,Z);
xlabel('x');
ylabel('y');
zlabel('f(x,y)');

title('The 3D Plot of f(x,y)');
end
% 定义适应函数
function z = fitness_func(x,y)
    z = 3*cos(x*y) + x + y^2;
end